# 设置工作目录  
setwd("F:\\r\\file")  

# 加载所需的R包  
library(dplyr)  
library(DESeq2)  
library(ggplot2)  
library(clusterProfiler)  
library(org.Hs.eg.db)  
library(enrichplot)  
library(gtable)  

# 安装缺失的包  
if (!requireNamespace("BiocManager", quietly = TRUE)) {  
  install.packages("BiocManager")  
}  
BiocManager::install("clusterProfiler")  
BiocManager::install("Biostrings")  

# 读取基因表达计数数据  
gene_expression_counts <- read.table("GSE42546_CleanedRawCounts.txt", header = TRUE)  

# 读取样本信息数据  
series <- read.table("GSE42546_series_matrix.txt", header = TRUE, fill = TRUE)  

# 数据预处理：设置行名为基因名，并移除多余的列  
rownames(gene_expression_counts) <- gene_expression_counts$GeneSymbol  
gene_expression_counts <- gene_expression_counts[, -1]  

# 数据CPM（每百万计数）处理  
data <- gene_expression_counts %>%  
  mutate(across(everything(), as.numeric)) %>%  
  mutate(across(everything(), ~ ./sum(.) * 1e6))  

# 对CPM数据进行对数变换并绘制直方图  
hist_data <- data %>%  
  as.data.frame() %>%  
  sapply(log2) %>%  
  as.data.frame()  

png('plot1_hist.png')  
hist(as.vector(unlist(hist_data)), breaks = 50, main = "Log2 Transformed CPM Distribution", xlab = "Log2(CPM)")  
dev.off()  

# 提取分组信息，这里假设分组信息在列名中  
group_info <- substr(colnames(data), 1, 1)  

# 构建样本信息表  
col_data <- data.frame(condition = factor(group_info))  
rownames(col_data) <- colnames(data)  

# 进行差异表达分析  
dds <- DESeqDataSetFromMatrix(countData = gene_expression_counts, colData = col_data, design = ~ condition)  
dds <- DESeq(dds)  
res <- results(dds)  
save(res, file = "res.RData")  

# 对比D组和C组，获取差异表达结果  
resDC <- results(dds, contrast = c("condition", "D", "C"))  
results_dfDC <- as.data.frame(resDC)  
head(results_dfDC)  
write.csv(results_dfDC, "deseq2_resultsDC.csv", row.names = TRUE)  

# 对比S组和C组，获取差异表达结果  
resSC <- results(dds, contrast = c("condition", "S", "C"))  
results_dfSC <- as.data.frame(resSC)  
head(results_dfSC)  
write.csv(results_dfSC, "deseq2_resultsSC.csv", row.names = TRUE)  

# 对比B组和C组，获取差异表达结果  
resBC <- results(dds, contrast = c("condition", "B", "C"))  
results_dfBC <- as.data.frame(resBC)  
head(results_dfBC)  
write.csv(results_dfBC, "deseq2_resultsBC.csv", row.names = TRUE)



# 添加颜色列，根据 log2FoldChange 和 padj 的值  
up_down_genes <- resDC  
up_down_genes$color <- "gray"  # 默认颜色为灰色  
up_down_genes$color[up_down_genes$padj < 0.05 & up_down_genes$log2FoldChange >= log2(1.2)] <- "red"  
up_down_genes$color[up_down_genes$padj < 0.05 & up_down_genes$log2FoldChange <= -log2(1.2)] <- "blue"  

# 保存上调和下调的基因  
up <- up_down_genes[up_down_genes$log2FoldChange >= log2(1.2) & up_down_genes$padj < 0.05, ]  
down <- up_down_genes[up_down_genes$log2FoldChange <= -log2(1.2) & up_down_genes$padj < 0.05, ]  
total <- rbind(up, down)  

# 绘图部分  
color_vals <- c("red", "gray", "blue")  
names(color_vals) <- c("red", "gray", "blue")  

ggplot(up_down_genes, aes(x = log2FoldChange, y = -log10(padj), color = color)) +  
  geom_point() +  
  theme_minimal() +  
  scale_color_manual(values = color_vals, guide = FALSE) +  
  labs(x = "log2 (fold change)", y = "-log10 (q-value)") +  
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey") +  
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), linetype = "dashed", color = "grey") +  
  theme(  
    axis.title = element_text(size = 16),  
    axis.text = element_text(size = 14)  
  )  

# 输出图片  
png('DC_volplot.png', width = 1000, height = 1000)  
print(last_plot())  
dev.off()

up_down_genes <- resBC  
up_down_genes$color <- "gray"  # 默认颜色为灰色  
up_down_genes$color[up_down_genes$padj < 0.05 & up_down_genes$log2FoldChange >= log2(1.2)] <- "red"  
up_down_genes$color[up_down_genes$padj < 0.05 & up_down_genes$log2FoldChange <= -log2(1.2)] <- "blue"  

# 保存上调和下调的基因  
up <- up_down_genes[up_down_genes$log2FoldChange >= log2(1.2) & up_down_genes$padj < 0.05, ]  
down <- up_down_genes[up_down_genes$log2FoldChange <= -log2(1.2) & up_down_genes$padj < 0.05, ]  
total <- rbind(up, down)  

# 绘图部分  
color_vals <- c("red", "gray", "blue")  
names(color_vals) <- c("red", "gray", "blue")  

ggplot(up_down_genes, aes(x = log2FoldChange, y = -log10(padj), color = color)) +  
  geom_point() +  
  theme_minimal() +  
  scale_color_manual(values = color_vals, guide = FALSE) +  
  labs(x = "log2 (fold change)", y = "-log10 (q-value)") +  
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey") +  
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), linetype = "dashed", color = "grey") +  
  theme(  
    axis.title = element_text(size = 16),  
    axis.text = element_text(size = 14)  
  )  

# 输出图片  
png('BC_volplot.png', width = 1000, height = 1000)  
print(last_plot())  
dev.off()

up_down_genes <- resSC  
up_down_genes$color <- "gray"  # 默认颜色为灰色  
up_down_genes$color[up_down_genes$padj < 0.05 & up_down_genes$log2FoldChange >= log2(1.2)] <- "red"  
up_down_genes$color[up_down_genes$padj < 0.05 & up_down_genes$log2FoldChange <= -log2(1.2)] <- "blue"  

# 保存上调和下调的基因  
up <- up_down_genes[up_down_genes$log2FoldChange >= log2(1.2) & up_down_genes$padj < 0.05, ]  
down <- up_down_genes[up_down_genes$log2FoldChange <= -log2(1.2) & up_down_genes$padj < 0.05, ]  
total <- rbind(up, down)  

# 绘图部分  
color_vals <- c("red", "gray", "blue")  
names(color_vals) <- c("red", "gray", "blue")  

ggplot(up_down_genes, aes(x = log2FoldChange, y = -log10(padj), color = color)) +  
  geom_point() +  
  theme_minimal() +  
  scale_color_manual(values = color_vals, guide = FALSE) +  
  labs(x = "log2 (fold change)", y = "-log10 (q-value)") +  
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey") +  
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), linetype = "dashed", color = "grey") +  
  theme(  
    axis.title = element_text(size = 16),  
    axis.text = element_text(size = 14)  
  )  

# 输出图片  
png('SC_volplot.png', width = 1000, height = 1000)  
print(last_plot())  
dev.off()

res<-as.data.frame(res)
ending1=res[abs(res$log2FoldChange)>log2(1.2) & res$pvalue<0.05,]
taggene=rownames(ending1)
#富集分析
gor=enrichGO(gene = taggene,
             OrgDb = org.Hs.eg.db,
             keyType = "SYMBOL",
             ont = "all",
             pAdjustMethod = "BH",
             pvalueCutoff = 0.05,
             qvalueCutoff = 0.05)

png(filename = "D:ot.png")
dotplot(gor)
dev.off()